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ABSTRACT 


There are currently no well-established numerical techniques for 
evaluating the bit error probability performance of a Satellite Communication 
System that includes: 

— Uplink and downlink noise 

— Uplink interference 

— Transponder AM/AM and AM/PM nonlinearities 

In this report we present new computational techniques that efficiently compute 
these bit error probabilities when only moments of the various interference ran- 
dom variables are available. The approach taken is a generalization of the well- 
known Gauss-Quadrature rules used for numerically evaluating single or multiple 
integrals. In what follows, we develop the basic algorithms, show some of its 
properties and generalizations, and describe its many potential applications. 

Some typical interference scenarios for which the results are particularly 
applicable include: 

— Intentional jamming 

— Adjacent and co-channel interferences 

— Radar pulses (RFI) 

Multipath 

— Intersymbol interference 

Wliile the examples presented stress evaluation of bit error probabilities in 
uncoded digital communication systems, the moment techniques can also be applied 
to the evaluation of other parameters, such as computational cutoff rate under 
both normal and mismatched receiver cases in coded systems. Another important 
application is the determination of the probability distributions of the output 
of a discrete-time dynamical system. This type of model occurs widely in 
control systems, queueing systems, and synchronization systems (e.g. discrete 
phase-locked loops) . 
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I. 


Introduction 


Our primary motivation for investigating the moment techniques presented 
here is the numerical evaluation of satellite communication system performance. 
These systems typically possess transponders which exhibit such nonlinearities as 
AM/AM and AM/PM conversion and are further corrupted by a combination cf uplink 

and downlink noise, and various interference signals such as those due to: 

— Intentional j amming 

— Ad j acent channels 

— Radar pulses 

— Multipath 

— Intersymbol interference 

— Co-channel interferers 

It is often difficult to have a complete statistical characterization of 
these interference signals. Some moments, however, are often easily computed 
based on some simple models of the various interference signals. Hence, given 
the available moments, we should desire a technique by V7hich one could achieve 
an approximate performance evaluation. 

The particular moment technique presented here is based on the solution 
to the classical ’^Hamberger Moment Problem" as discussed in Krein (Ref. 1). This 
solution has previously been applied to linear communication channels by 
Benedetto, Biglieri, and Castellan! (Ref. 2) and Yao and Biglieri (Ref. 3). It is 
also knoxm to be a generalization of the well-known "Gauss-Quadrature Rules" for 
numerically evaluating integrals (Ref. 4). We present here some new algorithms 
for solving the basic moment problem and then generalize them to complex and 
multi-dimensional random variables. Although our primary application is motivated 
by the evaluation of satellite communication system performance, there are numer- 
ous other practical applications of this moment technique which shall be discussed 
at the conclusion of this report. 

In Section II, we examine the transponder satellite model and motivate the 
need for developing a computationally efficient numerical technique for evaluating 
the bit error probability of such systems. The moment technique will have appli- 
cations to all types of signal modulations including the new bandwidth efficient 
modulations such as MSK, SQPSK, CPFSK, and TFM (Refs. 5-7). Section III dis- 
cusses the basic assumptions and statement of the classical one variable moment 
problem. Section IV presents a solution to the moment problem using the 




1 


Berlekamp-Massey algorithm (Ref, 8) and an accompanying root-finding algorithm 
along with some numerical examples illustrating their use. Section V presents an 
efficient algorithm for computing , the moments of a sum of independent random 
variables in terms of their individual moments. Section VI presents some basic 
existence theorems concerning the solutions to the moment problem. Generaliza- 
tions to complex random variables and pairs of correlated random variables are 
given in Section VII. Section VTTI shows how to solve the moment problem given 
some constraints on mass points. The accuracy of this moment technique as well 
as the derivation of bounds on the approximation error are presented in Sec- 
tion IX. Finally, various applications are discussed in Section X along with 
conclusions. 


Another algorithm which can be applied to this problem is the Euclid algorithm 
whose relation to the Berlekamp-Massey algorithm is discussed in Ref. 9. 
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II. Transponder Satellite Channel 


Typically, a transponder satellite is modelled as shown in Figure 1 where 
we define 

x(t) = transmitted uplink signal 
i(t) = uplink interference signal 
n^(t) = uplink noise 

r(t) = x(t) + i(t) + ~ signal entering the satellite system 

BPF = bandpass filter 
a(t) = signal entering the TWT 
TWT = traveling wave tube amplifier 
ZF = zonal filter 
z(t) = satellite downlink signal 
n,(t) = downlink noise 

Q 

y(t) = signal received at the ground station (2.1) 

We now illustrate the problem of performance evaluation for the above satel- 
lite channel when the modulation is coherent binary phase shift keying (BPSK) for 
which 


x(t) 


s(t) , "0" data bit is sent 

-s(t) , "1" data bit is sent 


( 2 . 2 ) 



Figure 1. Satellite Channel 
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where 


s(t) = /2P cos oj t ; 0 < t < T 

0 


P = transmitter power 


(2-3) 


We assume the BPF is ideal in that It limits the satellite input signal 
r(t) to the signal space generated by the pair of quadrature basis functions 


Hence 


where 






os ai^t 


= - /§ sin 0 ) t r 0 < t < T 
s -vj i 0 


(2.4) 


a(t) “ r (j) (t) -h r (p (t) 
c c, s s 


(2,5) 


r = 
c 


r(t)(j) (t) dt 
c 


X + i + n 
c c uc 


r = / r(t)(j) (t) dt = X -f i + n 

s f ^ s s us 


( 2 . 6 ) 


are the projections of r(t) on these basis coordinates. 
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Here for BPSK we have x = 0 and 

s 


X 

c 



"0" data bit 
”1” data bit 


(2.7) 


while 

i^, ig are the quadrature components of the interference signal 

%s independent components of the uplink additive white 

Gaussian noise. 

The bandpass filter’s function is to filter out all noise and interference 
outside this signal space (spectrum) without distorting the signal. We have 
assumed the bandpass filter works ideally. 

Next we define the envelope 


R ^ 



and phase 





of the signal a(t); i.e., 


a(t) = r^(j)^(t) + rg(j)g(t) 


= R cos + n] 


( 2 . 8 ) 


(2.9) 


( 2 . 10 ) 
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The TWT is assumed to create AM/AM and AH/PM nonlinear conversions which 
are mathematically described by the zevv menujj' functions of the input envelope R; 
viz. , 

f (R) = AM/Al'f nonlinearity 
g(R) = AM/PM nonlinearity 

Thus, the TWT output follov/ed by a zonal filter is given by 

z(t) = f(R) cos [(Opt + g(R) + n] 


= aJI' [g(R) + n] 'f'j.Ct) 


+JI f(R) sin [g(R) + n] <J>g(t) (2.11) 

In general, we ai' 5 sume a conventional/ground station receiver based on the 
ideal additive white Gaussian noise channel. With few exceptions, it is usually 
impractical to design special receivers for each channel. The conventional 
receiver is modelled as in B'’igure 2. 



Here we have 


Figure 2. Ground Station Receiver 




( 2 . 12 ) 
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( 2 . 13 ) 


+ n - g] 


g = receiver phase reference 


and the decision or demodulation rule 


decide ”0" if y >0 

decide ”1” if y ^ 0. 

•'c 


( 2 , 14 ) 


Suppose the data bit is sent |_5c(t) = s(t) = 

A PT is the energy per symbol. Then, given the conditional error probabil- 
ity is 


P„ (z ) “ Prob 




X = /e^ 
C 



( 2 . 15 ) 


where is the single-sided noise spectral density of 


the downlink noise n^(t). 


The average error probability is then 


E. 


= E Ip 

j E, 




( 2 . 16 ) 


We assume a phase-locked loop tracks the long time a/erage phase of the satellite 
output signal. 

The function 


Q(x) = 


V 2 7T 


exp(-y )dy 


is the well-known Gaussian probability integral. 
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whuro E {•} denotes the expi'ctution over the probability density function of z^. 

Since from (2.13) toge,ther with (2.8) and (2.9). z is a function of r and r 

c c s 

with now 


r = + i + n 

c s c uc 


r = i + n 
s s us 


then, equivalently 


/ 9 I 7 > 

1 ir * ^ * J 


t = F / /e~ + i+u ,i+n \ 
C \ ® uc S US J 


(2.18) 


for some known function F (.,.). 

Hence, the average bit error probability has the form 


\ ^ r^o 


F ( /e~ + i + n , i + n 

\ S C uc S US 


n \ 


(2.19) 


where E {•} is now the expectation over the random variables 1 ,i , n , and n . 

C S Lie.. Llo 

Another form for this error probability can be had by first defining the 
complex random variable 


where 


W = /t 


/E + i + 
s c 


n ^ + j (i^ + n ) 
uc I s us 


= Re 


in 


j - /“I 

R = jw| 


( 2 . 20 ) 


( 2 . 21 ) 
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Then we have for some known function G(*) the error probability 



E jG(W)j 


( 2 . 22 ) 


The key to evaluating the bit error probability for 
technique with the general satellite channel involves the 


the BPSK modulation 
evaluation of 


P 


I "c 


[F(r 




■= E[G(W)] (2.23) 

wheivi r^ and r^ are, in general, two correlated real random va:.lables. This 

requires knowledge of the joint probability distribution p(r ,r ) of the pair of 

c s 

random variables r^ and r^. This is not always available and even when we have 
it, we still have to evaluate the double integral 

In practice, it is often easier to obtain'some joint moments 





?.,m = 0,1,2, ••• , N 


'r complex moments 


= E(W“) ; m = 0,1,2,--*, N 


(2.25) 


(2.26) 
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The remainder of this report examines how we can use available moments and 
obtain an approximation to E E[G(W)], In particular, we 

describe a way of obtaining discrete approximate probability distributions of the 
form* 


Pr(W = ^ “ 1,2,' • • , V 


based on morents of the complex random variable W as in (2.26) or 

A 

“ Pjl 5 A = 1,2,*-*,V 


(2.27) 


(2.28) 


based on joint moments of two real random variables r^, r^ as in (2.25). These 
approximate distributions then yield approximations to (2.23) in the form, 


E 





(2.29) 


and 


E 



E 

0=r1 




(2.30) 


Equations (2.29) and (2.30) represent generalizations of the Gauss-Quadrature 
technique which is often applied to numerically evaluate double Integrals of the 
form in (2.24) when p(r^,r^) happens to be of a Gaussian nature. 

Although we limited this example to BPSK, the approach outlined here 
applies equally well to coherent MPSK for all integer M and also to other general 
bandwidth efficient modulation techniques. 


*The hat is used to denote the word "approximate." 
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The Classical One Variable Moment Problem 




Let X be a random variable (continuous or discrete) and suppose we only 
know its N + 1 moments 


= E(X^); k = 0,1,2,‘ • • , N O.l) 

where y^ = 1. We want to find an approximation to the true probability distri- 
bution of X in the form of a discrete probability distribution. The classical 
moment problem is to fln'l the smallest number of points x^, x^, and 

weights so that the approximating distribution 

Pr|x = x^l = 0)^ ; £ = l,2 ,---,v (3.2) 

satisfies the given moment constraints, 

V 

= E(XjP = lo^x^ ; k = 0,1,2,- • • , N (3.3) 

£=1 

Suppose we start by assuming that Pr(X = x^) = Z = 1,2,---, v is the 
true distribution so that 


“k ■ ■ Z) 

£=1 


k 


(3.4) 


is true for all values of k = 0,1, *••• Next define the polynomial 


V 

C(D) = (1 - Dx^) 

)l=l 


= c. 


CiD 


^2^ 


+ c D 
V 


(3.5) 
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where Cq = 1. Also consider, for arbitrary n, the relation 


vr^ 

j=0 j=0 \ £=1 


\j=0 


X 


-J 


■ I] 

£=1 


(3.6) 


since is a root of C(D). Then recalling the fact that c^ =1, (3.6) can be 
written in the alternate form 



(3.7) 


This form of the relationship between moments allows us to interpret moments as 
outputs of a real field linear feedback shift register as shown in Figure 3. 



Figure 3. Moment Generating Linear Feedback 
Shift Register 
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Note that, although at this point, we do not know the points x^, x^, *•*, 
nor the polynomial C(D) given in (3.5), we have the interpretation that the 
given N + 1 moments of (3.1) are generated by some linear feedback shift 
register with feedback coefficients that specify this polynomial. This is a new 
interpretation or formulation of the classical moment problem. 

Next define the polynomial 


P(D) 


V V 


■£-.n 

Z=1 j=l 




(1 - 


= Pg + PjD + p^O^ + • • • + p^_^ 

Then the moment generating function polynomial 


P(D) A 


00 


k=0 


D 


k 



D 


k 



(3.8) 


(3.9) 
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when multiplied by the polynomial C(D) yields the relation 


> 


V 

p(D)C(D) =^0)^ 


V 


O'" 

•1=1 




DXj) 


(3.10) 


= P(D) 


By equating terms with equal powers of D, the coefficients of P(D) are given as 
follows: 


Po ■ “o 


Pi ■ + °1 Pq 


^2 ^2 ^1 ^1 ' ^2 ^0 


v-1 


P . + Ct P „ + 

V-1 1 v-2 


+ Pv-l^O 


(3-11) 


Thus, given the polynomial C(D) and the known moments of X, we can easily obtain 

the polynomial P(D). Given these two polynomials we show next how the weights 

o)„, •••, 0 ) are easily found. 

X ^ V 
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Assume we have polynomials C(D), P(D) and the reciprocals of the roots of 
C(D) which are the points x^. Then, from (3.5) 


C»(D) C(D) 


Note that 


and from (3.8) 


Thus, 


V V 


j=i 


(1 - 


DXj) 


c, + 2c„D + 3c„D^ + • • • + VC 
1. 2 3 V 


(3.12) 


c 


- X, 


=1 


K ^- 


"kS) 


(3.13) 


P(x, ^ 1 = 0 ). 


fl (' ■ 


j=l 


(3.14) 



-ll 

J 


C'(k 

1 — 1 
I 

, iC X,Z, ,v 


(3.15) 
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From the above relationships, we see in summary that the classical moment 
problem is solved by first finding the shortest length linear feedback shift 
register that generates the given N + 1 moments. This feedback shift register 
is specified by the polynomial C(D) whose roots have reciprocals which are the 
desired probability mass location points X 2 » •'*, x^. Next obtain P(D) from 
(3.11) and the probability mass values **'» from (3.15). 

In the next section, we describe two basic algorithms, namely the 
Rerlekamp-Massey algorithm which enables one to find the polynomial C(D) and an 
algorithm to find the roots of C(D). 
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IV . The Be.rlekamp-Massey Algorithm (Ref. 8) 


Given y^, y^, ••• , y^ the Berlekamp-Massey linear feedback shift register 
synthesis algorithm is a technique for finding a smallest length feedback shift 
register that generates y^, y^, •••, y^ and is described by the polynomial 

v> 

C(D) = 

1=1 


•■= Cq + + • • • + 


(4.1) 


The following is a step-by-step description of this algorithm. 
Define the following variables: 

(a) m. n, £ integers 

(b) b, d real numbers 


(c) C(D),B(D), T(D) 


C(D) = 1 -f Cj^D + c^D + ••• + c^D 


polynomials in D 
£ 


Step 1: Input moments 

Pq = 1. ^2’ ” ’ ’ 

Step 2: Set initial conditions 

C(D) = 1, B(D) = 0, T(D) = 1 

m=l, n = 0, )l = 0, b = l 

Step 3: Compute 

d = y + c,y - + c„y „+•••+ c y 

n 1 n~l 2 n-2 Z n-x. 
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Step 4: 


If d ~ 0, then 
m + 1 m 
and go to Step 7. 

Step 5: If d 0 and 2£ > n, then: 

C(D) - ^ d"‘b(D) ->• C(D) 

m + 1 ->■ m 
and go to Step 7. 

Step 6: If d ^ 0 and 2H < n, then 

C(D) ->• T(D) 

C(D) - f D™B(D) C(D) 
b 

n + 1 -<!.•+ £ 

T(D) -5- B(D) 
d b 
1 m 

Step 7: n + 1 ->• n 

Step 8: If n = N + 1, stop. Otherwise go to Step 3. This algorithm results 

in C(D) of (4.1) 

and 

P(D) = Pq + P^D + p^D^ + ... + (4.2) 


The notation "A-^B" means replace B with A. 
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where from (3.11) we have 



and because of the distinct condition, 


X . I only if x , = -x . . 
3' 1 3 


In general, the linear feedback shift register relationship 



k = 2v + 1, 2v + 2, • • • 


(4.5) 


for some initial conditions is satisfied by outputs of the 

form 



k = 2v + 1, 2v + 2, “ • 


(4.6) 


Section VI will prove these reciprocal roots are distinct and real. 
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with arbitrary coaff iaients a^, 'it 2 , •••, To so« this, substitute' (4.6) into 

(4.5) which prcjduces 


V 

■ 2-J 


,i=l 




E k /V' -r 

a.x. I > c.x."' 

1 1 \ jLu J 1 i 
j=l 


(4.7) 


i=i 


Since from (4,1) the quantity in parentheses can be identified as C - Cq 

which from the factored form of (4.1) is seen to have value --I, then (4.7) 
immediately reduces to (4.6). 

Note that only when the initial conditions of the feedback register are 
set to the given moments of the random variable X do v/e necessarily have 
ai~ai^; i = 1, 2, ' • ’ , v. Here we consider arbitrary initial conditions for the 
linear feedback shift register defined by C(D). 

Kc-v-t consider for some k > 2v, 


V 

E k 

a.x, 
1 1 


i=l 


k ; “2 / 2 

= a.x. <1 + — I — 
1 1 i “1 \^1 


k 


+ 


“1 


+ . . . + 


a /X 
V / __y_ 

“1 Vi, 


(4.8) 
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Def iitc 



and note that if the magnitudes of the x^'s are ordered as in (4.4), then 

lim F(k) = 1 (4.10) 


where convergence is primarily determined by the term, 



Here we can take the ratio of the feedback shift register outputs 


^k+l ^ F(k + 1) 
""l F(k) 


(4.11) 


and find the first reciprocal root by 


X, 


lim 

k->“ 



(4.12) 


If we have |x^| 
the form 



V then x^ = -x^ and F(k) has 


“2 k “3 

F(k) = 1 + / (-1) + — 

“1 “1 


fx. 


+ 



(4.13) 
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which oscillate.B between the. limits 
have 


“2 

1 *: as 

”1 


k increases to infinity. . Thus, we 


lim 

k-x» 


F (k + 2 ) 
*F(k) 


(4.14) 


and the ratio 


^‘k+2 ^ 2 F( k + 2 ) 

F(k) 


(4.15) 


yields the first two reciprocal roots 



(4.16) 


Note that we need to choose the initial condition of the feedback shift 
register such that 0 for k > 2v. 


An 


efficient way to generate is to define 


Then 



(4.17) 


r-x 


A A , 
r r-1 


r--i+l 


(4.18) 
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Now rGcalllng from (3.7) with k = n that 


-'^k ■ "2\-2 + ■ ■ ■ %"k-v 


(4.19) 


then dividing by ^ we get 


c- + - + c„A A, ^ 4- 

A- 1 2 k*"! 3 k“l k-2 

k 


4- c A, - A- ^ 
V k“l k~2 




This recursion relationship together with (4.12) and (4.17) gives 


lim ^ 
k-x» k 


(4.21) 


provided |x^| > 2,3, v. Alternately using (4.16), the first two 

reciprocal roots become 



(4.22) 


The procedure for finding the remaining reciprocal roots is outlined as 
follows. ' Suppose we find as described above for the case where |x^| > 

I = 2,3, V. Then, we remove the corresponding factor from C(D) and define 

a new polynomial 
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V 

C^^^(D) = J~J (1 - Dxp 

il=2 


C(D) 

1 - Dxj_ 


- .(1) 


= c 


+ + 


v-1 


(4.23) 


From the relation 


C(D) = (1 - (D) 


(4.24) 


we equate the coefficients of equal powers of D and obtain the relations 




( 1 ) 


( 1 ) ( 1 ) 

^1 *^1 " ^ 1^0 


( 1 ) ( 1 ) 

^^2 ^2 “ ^ 1^1 


( 1 ) ( 1 ) 
c _=c '-XtC ' 
v-1 v-1 1 v-2 


( 1 ) 

C = - X, C , - 

•; 1 v-1 


(4.25) 
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•'U7 f. 


or equivalently 


^0 0 ^ 


= c + X 
1 1 10 


( 1 ) 




(1) . (1) V 

c i=c - +X-C ;!= 

v-1 v-1 1 V“2 


(4.26) 


.( 1 ), 


The recursive relations in (4.26) define the polynomial (D) . 

If we have a pair of reciprocal roots such that then we first 

remove both of the corresponding factors from C'(D) and define a new polynomial 


V > 

= j~J (1 - Dxp) 
ii=3 


COl) 

(1 - Dx^) (1 - DX 2 ) 


C(D) 


1 - 


2n2 

XiD 


= ^ cf )0 . 




(4.27) 


25 



From the relation 


we obtain 


C(D) = 





,( 2 ) 

'0 

,( 2 ) 

'1 

,( 2 ) 


= c. 


. 2 ( 2 ) 


(4.28) 


( 2 ) . 2 ( 2 ) 

v-2 v-2 1 v-4 


( 2 ) 

The recursive relations in (4.29) define the polynomial C (D) 


(4.29) 


The above approach for finding the largest magnitude reciprocal root or 
roots is then applied to the new polynomial C^^^(D) or C^^^(D) to find the next 
largest magnitude reciprocal root or roots. This procedure is continued until 
all reciprocal roots have been found. 


The following is a summary of the reciprocal root-finding algorithm just 
described. 


Root Finding Algorithm 

Assume moments y^, •••, y^ are used in the Berklekamp-Massey algorithm 

to find the polynomial 


C(D) == Cq + c^D + c^D 


4 - 


c D 

V 
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where - 1, The roots of C(D) are unique and real. The following algorithm 
finds their reciprocals; 


Step 1: Input 


C- y , • • • c ; N and e 
1 ^ V e 


Step 2: Set 


Step 3: 


Step 4: 


^2 ^ v-1 “ ^ 


Set 

k = V 
Compute 


z, = c, + , + c^X, .X, ^ 1- c. X,_ ^X. 


^k 1 2 k-1 3 k-1 k-2 


Step 5: Compute 


^k z, 
k 


Step 6; Compute 


V k~l k-2 k-v+1 


■ Vk-i 


Step 

7: 

If l\ - \_il 

+ 

Step 

8: 

If k > N^, go 

to 

Step 

9: 

k k + 1 and 

go 

Step 

10: 

i£ hk - vJ 

> 

Step 

11: 

Set 



X = -z, 
V k 


and e are convergence parameters. 


27 


Step 12: 


Step 13! 

Step 14: 

Step 15: 
Step 16: 


Step 17: 


c- + X c. 

1 V 0 1 


■=2 * =2 


C , + X C o c , 

v-1 V v-2 v-1 


V v-1 
If V = 1 , 

set Xj^ = -c^ and stop. 
Go to Step 2. 


X = /tT 

V k 


X . = -/tT 

v-1 k 


=2 + * '2 


“=3 + * =3 


Step 18: 

Step 19: 
Step 20: 

Step 21: 
Step 22: 


c „+xc ,^c „ 

v-2 V v-4 v-2 


V V - 2 

If V = 1, set x^ = -Cj and stop. 

If V = 0, stop. 

Go to Step 2. 

Declare ill-conditioned and redo Berlekamp-Massey algorithm with 
moments y^, •••, 
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a 


If the original moments are in fact not true moments, then 

the Berlekamp-Massey algorithm can result in a polynomial C(D) whose roots are 
complex. This can be caused by various errors in computing these moments, as 
well as possible roundoff errors in the above algorithms. Step 22 attempts to 
detect such problems. 

Typically small changes in the coefficients c^, c^, •••, can cause large 
changes in the roots of C(D), particularly the larger roots. The smaller roots 
of C(D) are generally more stable. This means that for the recipi^ocal roots 
Xf, x^, •••, x^, the larg'er magnitude points tend to be more stable. 

To reduce roundoff errors in the Berlekamp-Massey algorithm, it helps to 
control the dynamic range of the moments 

= E(X^) ; k = 0,1,2, ••*, N (4.30) 


by defining 


Y = pX 


(4.31) 


with moments 


Pk'P) 5 ^ “ 0,1,2, ••*, N (4.32) 

If we apply the Berlekamp-Massey and the root-finding algorithms to the moments 
of Y = pX, then the resulting mass location points y^, y^ , •••, y^ are related 
to the desired points x^, X 2 , by 


X 


I 



I = 1,2, • • • , V 


(4.33) 


The weights •••> remain the same in both cases. Here p can be 

selected to control the dynamic range of the input moments to the Berlekamp- 
Massey algorithm. A good choice is governed by the condition 

y2(P) = 1 (4.34) 
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or 



(4.35) 


We conclude this section with a 1 numerical examples to illustrate the 
use of the algorithms just discussed. The examples chosen will correspond to 
probability distributions for which all moments are known. Thus the end products 
of applying the foregoing algorithms will serve as verification of well-known 
Gauss-Quadrature results for these distributions (Ref. 4). 

As a first example, consider a zero-mean Gaussian probability density 
function for which 


= l-3'5 ••• (2n-l) A (2n - 1) I : 
zn 

li 2 ^_l = 0 ; n = 1,2. ... 


(4.36) 


Assume for the purpose of this example that only the first ten moments in 
(4.36) are known. Then, using these as an input, the Berlekamp-Massey algorithm 
proceeds step-by-step, as follows: 

Step 1: (N = 9) 

Uq = 1, = 0, == 1, Pg = 0, p^ = 3, 

p^ = 0, Pg = 15, p^ = 0, Pg = 105, Pg = 0 

Step 2; 

C(D) = 1, B(D) = 0, T(D) = 1 
m=l, n=0, £ = 0, b = l 

Step 3: 
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Step 6: 


Step 7: 

Step 8: 
Step 3: 

Step 4: 

Step 7: 

Step 8; 
Step 3: 

Step 6: 




(2i = n) 

T(D) = 1 

C(D) =1=1- (0)D ; = 0 

1 = 1 
B(D) = 1 
b = 1 
m = 1 

n = 1 
(n < 10) 


m = 2 

n = 2 
(n < 10) 

d = P2 + = 1 

(2£ = n) 

T(D) = 1 

C(D) = 1 - = 0, C2 = -1 

5, = 2 + l- l= 2 
B(D) = 1 
b = 1 
m = 1 


31 


Step 7: 


Step 8 : 
Step 3: 


Step 4: 


Step 7: 


Step 8 : 
Step 3: 


Step 6 : 


Step 7: 


Step 8 : 
Step 3: 


n = 3 
(n < 10 ) 

d = M 3 + c^M 2 + c^M^ = 0 
m = 2 


n = 4 
(n < 10 ) 


,=0 

<1 - ^4 + ■:/'3 + 


(- 1 ) ( 1 ) 


( 2 £ = n) 


T(D) = 1 - 

C(D) = 1 - - 2D^ = 1 - 3D^; = 0, = -3 

£ = 4 + l- 2 = 3 
B(D) = 1 - 


b = 2 
m = 2 


n = 5 
(n < 10 ) 


d = Me 


+ ^/^4 4^3 ^ 4^2 = ° 
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Step 4: 


Step 7: 


Step 8: 
Step 3: 


Step 6; 


Step 7: 


Step 8: 
Step 3: 


Step 4: 


Step 7: 


m = 2 

n = 6 
(n < 10) 

d - ^ ‘^/^3 ^ ^ 

15 (-3) (3) 

(2Z = n) 

T(D) = l‘- 3D^ 

C(D) = 1 - 3D^ - (1 - D^) 

2 4 

= 1 - 6D 4- 3D ; = 0, = 6, - 0, = 3 

£ = 6 + 1- 3^4 

B(D) = 1 - 3D^ 

b = 6 

m = 1 

n = 7 
(n < 10) 



m = 2 


n = 8 
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tJ 


i 



Step 8: 
Step 3: 


(n < 10) 




Step '6: 


Step 7: 


Step 8: 
Step 3: 


Step 4: 


Step 7: 


Step 8: 


d 


^8 

ToT 



‘’3/5 ■*' 

(-6)(15) 


(3H^ 


24 


(2A = n) 

T(D) = 1 - 6D^ + 3D^ 

C(D) = 1 - 6D^ + 3D^ - (1 - 3D^) 

= 1 - lOD^ + ISd"^; = 0, = -10, = 0, c,^ = 15, = 0 

£--=8 + l - 4 = 5 
B(D) = 1 - 6D^ + 3D^ 
b = 24 


m = 1 


n = 9 
(n < 10) 



m = 2 


n = 10 

(n = 10) . Stop. 
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The resHilting linear feedback shift registe^r analogoUvS to Figure 3 is 
illustrated below: 



Figure 4. Moment Generating Linear Feedback Shift Register for 
Ten Gaussian Moments 


The corresponding generating polynomial is 

C(D) = 1 - lOD^ + 15D^ 


(4.37) 


which is the desired result. 

Note that the last value of S, [the order of the polynomial C(D)] computed 
by the algorithm is H - 5. Thus, since (4.37) is only a fourth order polynomial 
in D, we immediately conclude that 


c^ = 0 (4.38) 

Equivalently, from the factored form of C(D) in (4.1), (4.38) tells us that 
one reciprocal root has value zero; i.e.. 



(4,39) 
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The remaining roots can easily be obtained by solving a quadratic equation or 
applying the root-finding algorithm. In the former case, let Z = D in (4.37) 
and equate the result to zero, namely. 


1 


lOZ + 15Z = 0 


(4.40) 


whose solutions are 


10 ±Ao 
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.544151844, .122514823 


(4.41) 


or 


D = ±.737666486, +.350021175 


(4.42) 


Finally, the corresponding reciprocal roots are 


3 = +1.35562618 

(4.43) 

X, . = ±2.856970014 
4,5 


Before showing how the root-finding algorithm can be used to approach the 
results in (4.43), we shall finish the solution for the approximating probability 
distribution by finding the five weights oi^. From the coefficients 

of C(D) and the given moments, (3.11) allows us to compute the coefficients of 
the polynomial P(D) which for this case becomes 



1 0 0 0 0 

0 1 0 0 0 

-10 0 10 0 


0 -10 0 1 0 


15 0 -10 0 1 


C 


1 

0 

1 

0 

3 


u 


(4.44) 
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or 


Pn “ Pi = Po = P-i = Pa “ ^ 


P(I)) = 1 - 9D^ + 


Differentiating (4.37) with respect to I) gives 


(4.43) 


C (D) = ~20D + 60D' 


(4.46) 


Finally applying (3.15), we get the distribution weights 


- . -ife'L 

K') 


2 4 

8 - 9x- + X. 
k k 

60 - 20x? 


(4.47) 


or, using (4.43) 


^2 3 = .222075922 


^ “ .011257411 


(4.48) 


Clearly, if we try to apply (4.47) to the reciprocal root = 0, we get the 
result 'll = -8/60 which is meaningless since probability distribution weights 

K. 

cannot be negative. Thus, whenever one of the reciprocal roots is zero, we must 
determine its corresponding weight from the usual normalization constraint on 
probability distributions, namely, 

V 


E 

k=l 


m, = 1 
k 


(4.49) 
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Letting V = 5 and subvStituting (4. 48) In (4.49) gives the remaining desired 
result, namely, 

= .533333333 (4.50) 

To check with the result given in Ref. 4 for Gaussian-Herraite Quadrature, 
we need to divide the reciprocal roots {x^} of (4.39) and (4.43) by /2 and multiply 
the weights {o)^} of (4.48) and (4.50) by t^. When this is done, we obtain exact 
agreement with the tabulations for n = 5 in Appendix B of Ref. 4 of page 343. 

We. now demonstrate how the root-finding algorithm can be used to rapidly 
approach the results found in (4.43) by solution of a -quadratic equation. 


Step 1: 


Step 2: 


Step 3: 


c^ = 0, C2 = -10, Cg = 0, c^ = 15 


^1 ^2 ^3 ^ 


k = 4 


Step 4: 


Step 5: 


Step 6: 


Step 9: 


z, = -10(1) + 15(1)(1)(1) = 5 




= -5 


k = 5 
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Step 4: 


Zj. = + 15(-i)(l)(D = -1 


Step 5: 


Step 6: 


Step 9: 



1 



k = 6 


Step 4: 

Zg = -10(1) +15(1) (1) = -13 

Step 5: 


X 


6 


13 


Step 6; 



Step 9: 


k = 7 


Step 4: 

.7 - -lo(i) + 15 (i) (1) (- 

Step 5: 


= 1 


Step 6: 


T 


7 



= 13 
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rH ( LT) 





Step 7: 


|Ty - Tfil + I'i’e - T^l = 0 + 18 = 18 


t ep 9: 


2g = -10(1) + 15(1)(^) (1) = 


115 

13 


X„ = 


13 


8 115 


^8 = 


{iis) 


( 1 ) 


115 

13 


ITg - I^l + |T, - T^l- 


115 


13 


-13 +0=7 


54 

13 


= 4.153846 


k = 9 




X,- 1 


T = 
9 


( 1 ) 


(iS) 


115 

13 


|Tg - Tgl + |Tg - T^| = 0 + 


115 

13 


- 13 


54 

13 


k = 10 


z^O 10(1) + 15(l)|^j^^ j(l) 23 

, =-23_ 

10 191 


191 


-■(3t)(l)-- 


jV 

Herein we avoid writing out the particular steps we are at since the sequence 
is always Step 4, Step 5, Step 6, Step 7, Step 9 until convergence is 
obtained. 
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"^ 10 - ^9 + ^ 9-^8 


1191 115 


,541806 


Notice how rapidly |Tj^ - + pk-l ~ "^k 2 converging. However, 

Zy^ - is not. Thus, ultimately the test in Step 7 will be satisfied, and 

so will the test in Step 10 which takes us to Step 16, namely the solutions for 
the two reciprocal roots of largest magnitude. Let us examine how close we are 
to the true results in (4.43) at this point in the root-finding algorithm. From 
Step 16, we have 

X. s 7 t77 = = 2.881726536 





= -2.881726536 


Comparing (4.51) with (4.43), we observe that after only 7 iterations of the 

algorithm, we are already quite close to the true result, namely 

X, . = ±2.856970014. 

4,5 

The next example chosen for Illustration is a uniform distribution; i.e.. 


p(x) = 


X < 1 


0 ; X > 1 


The moments of this distribution are easily found to be 


1 i K 

“k ■ 2 - ■ 


k + 1 


k odd 


; k even 
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Again let's start by assuming knowledge of only seven moments. Than, the 
Berlekamp-Massey algorithm proceeds as follows: 


Step 1: 


Step 2: 


Step 3: 


Step 6: 


Step 7 : 


Step 8: 
Step 3: 


Step 4: 


(N = 6) 

Pq = Pj^ = 0, s ’’ ^^3 J ’ *^5 ” “ y 

C(D) = 1, B‘(D) = 0, T(D) = 1 
m = l, n = 0, X. = 0, b = l 

d = Pq 1 
(21 = n) 

T(D) = 1 

C(D) =1=1- (0)D ; = 0 

1 = 1 
B(D) = 1 
b = 1 
m = 1 

n = 1 
(n < 7) 



m = 2 
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Step 7: 


Step 8: 
Step 3: 


Step 6: 


* 

Step 7: 


n = 2 
(n < 7) 


(25, = n) 

T(D) = 1 

C(D) = 1 - = 0, 

«, = 2 + l- l = 2 
B(D) = 1 



m = 1 


d = ^3 + ° 

m = 2 


1 

3 


n = 4 


d = 
T(D) = 

C(D) = 
Si = 


\ ^ l '^3 ‘^ 2*^2 


1 2 
1 - jD^ 


Jl 

(ir 

4 + 1 - 2 = 3 


1 - -m - 


1 

= i + 0 + 



45 



3 

5 


Here again we shall omit the step 


numbers until we reach n = 


7. 
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n = 5 


n = 6 


B(D) = 1 = 
b 

m = 1 

d = 1J5 + + C2U2 + c^y2 = 0 

m = 2 


d = + C3«3 . i + 0 + + 0 

= ^ 

175 



£=6 + l- 3 = 4 
B(D) = 1 - 



m = 1 

n = 7. Stop. 


Since the last value of I (namely £ = 4) in this case agrees with the order 

of the final polynomial C(D), there is no reciprocal root which has value zero. 

2 

The four reciprocal roots can be obtained as before by substituting Z = D in C(D) 
and solving the resulting quadratic equation. In particular. 


1 




= 0 


(4.54) 
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whose solutions are 


Z = 8.651483715, 1.348516283 


or 


D = +2.941340462, +1.161256338 


Finally, the corresponding reciprocal roots are 


2 1.339981044 

X. , = ±.861136312 
3,4 

Again the x>^eights of the approximating probability distribution are 
by substituting the given moments and the coefficients of C(D) in (3.11). 


or 


1 

o 


1 

0 

0 

0 

1 

Pi 


0 

1 

0 

0 

0 

P2 


6 

“ 7 

0 

1 

0 

1 

3 

j3_ 

f 

0 

6 

" 7 

0 

1 

0 

Po 

= 1, 

Pi = 

0, P2 = 

11 

21’ 

P 3 = 

0 


P(D) = 1 - 


Differentiating C(D) with respect to D gives 


C'(D) 


12 n . .12 3 
7 ° ^ 35 ° 


(4.55) 


(4.56) 


(4.57) 

found 
Thus , 


(4.58) 


(4.59) 


(4 . 60) 
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Finally applying (3.15), we get the distribution weights 



4 11 2 

\ " 21 
12 12 2 

35 " 7 


or using (4.57), these evaluate to 


2 = .326072569 

to- , = .173927423 
3,4 


(4. 61) 


(4.62) 


To check with results given in Ref. 4 for Gauss-Quadrature with constant 
weight function, we merely need to multiply the weights of (4.62) by 2. When 
this is done, we obtain exact agreement with the tabulations for n = 4 in 
Appendix A of Ref. 4 on page 337. 
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V. 


Computing Moments of S ums 


In many applications, we wish to compute the moments of the sum of 
independent random variables. An efficient algorithm for doing this when given 
the moments of the individual terms in the sum is presented here. This approach 
is due to T. C. Huang (Ref. 10). 

We assume the random variable X with moments 


= E(X*^) ; k = 1, 2, ••• 
has a moment generating function 


$(0)) = E(e“^) 


Using the expansion 


00 



k=l 


(5.1) 


(5.2) 


(5.3) 


the moment generating function is given in terms of the moments by 


CO 


$(o)) = 1 + 


E 



Next use the expansion 


in (1 + a) 



(-1) 


j+i 


(5.4) 


(5.5) 


47 


to obtain the form 




Jin $(w) = Jin 




k 



00 


E 

0 = 1 



X 


I 


(5-6) 


Here are the so-called semi-invariants of X and can be expressed as 

a weighted sum of the moments. 

We now determine an algorithm for computing the semi-invariants from the 
moments and vice-versa. 

Define 


Then , 


E(o)) 


0 = 1 


Jl 

TT^o 

Jl! Z 


(5.7) 


$(o)) 


^E(„) 


(5.8) 


has derivatives 


(w) 


■EM' 




k=0 


n 


j-1 „ . - 


n = 1,2, 


(5.9) 
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Since from (5.4) 


and from (5.7) 


♦ (0) . m 


E^-’''(0) = X. 


then evaluating (5.9) at u) = 0 gives us the desired relationship, namely* 


n ^yMm .X. 


'n^E(S) 


m .X . 

1 


or equivalently 




.m 

2 n-3 


Here (5.12) and (5.13) together with the initial condition 


™i ^1 


allows us to easily compute semi- invariants from moments and moments from 
semi- invariants , 


*Note: m^ = $(0) = 1. 
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Suppose now we have a sum of independent random variables 


f • • 


9 


X, 




i. e. , 


Y = + X 2 + • •• + Xj^ (5.15) 

and we wish to find the moments of Y defined by 

= E(Y^) ; k = 0,1,2, •••,N (5.16) 

when given the moments of the individual namely, 

niik = E(xJ) ; k = 0,1,2, N 

i = 1,2, • • • , L (5,17) 


We begin by defining a recursion equation analogous to (5,13) which relates 
the moments and semi-invariants of each random variable X^, namely, 


Xj - m. 
in in 


-E 

j=l 


m . 


:L,n-j 


n = 1,2, 


N 


(5.18) 


where 


X.. = m. - for i=l,2, 
il ii 

Next, recall from (5.7) and (5.8) that 


(5,19) 



I 

EVe I = .exp 


(5.20) 



Now consider the moment generating function of Y as follows: 


$y(o)) = E(e“^) 


' oiCXj^+X^H- 
e 




i=l 



(5.21) 


Thus, the moments of Y are obtained from a recursion relation Identical to (5.12) 
i.e, , 



(5.22) 
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where 


L 

AV 

1=1 


i£ 



f^,23) 


Figure 5 is a flow chart representation 



Figure 5. Moments of Sums 


which shows how the moments of Y are easily obtained from the moments of 
•••, X^. The procedure involves L transformations, T, from moments to semi- 
invariants using (5.18), taking the sum of these semi- invariants to obtain the 

-1 

semi- invariants of Y, and finally inverting the transformation T once, using 
(5.22) to obtain the desired moments of Y. 


As an example of the application of the results in this section, consider 
the important problem of assessing the performance of the satellite communication 
system modeled in Section II in the presence of multiple pulsed RFI sources. For 
the purpose of this example, we assume that each RFI source emits pulses with 
Poisson arrival times and the sources are independent of one another. Thus, for 
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til 

the i source, i = 1, 2, •••,!, the probability that n pulses occur in an 
Interval T is described by the distribution 


p(n) = e 





n = 0,1,2, •• • 


(5.24) 


where the mean of the distribution, is typically linearly related to T, i.e», 

y. = a.T (5.25) 

' 1 X 


We wish to characterize the moments of the discrete random variable corresponding 
to the total number of pulses in an interval T contributed by the L sources. 

The random variable corresponds to the number of pulses which arrive 
from source i in the interval T. Using the Poisson distribution of (5.24), we 
compute the moment generating function of as 



n=0 


0) / tAl T \ 

. y .e y , (e -1) 

^ X ‘ X X 

= 8 e = e 


(5.26) 


Using (5.6), we can immediately identify the semi- invariants of X^ as follows: 


^n $ (w) = y.(e ~1) 

X. ' X 

X 


00 


n=l 


00 



n=l 


(5.27) 
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or 


A. = Y, for all n 
in 1 


(5.28) 


Thus, for a Poisson process, we see that all the semi-invariants are equal to the 
mean of the process. 

Letting Y of (5.15) now correspond to the random variable characterizing 
the total number of pulses in the interval T contributed by the L sources, then 
we can immediately apply (5.22) and (5.23) to obtain its moments. Thus, 


X- - 


n 



for all n 


(5.29) 


and 



(5.30) 
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VI. 


Existence and Uniqueness of Solutions 

Given the moments y^, • • • , y^ of the random variable X, the Berlekamp- 
Massey algorithm finds the smallest number v and coefficients Cj^, c^, •••, 
such that 


y 


n 


-Z 

0=1 




n=v, v+1, N 


( 6 . 1 ) 


where if N is odd then* 


, N + 1 

V < — = — 


( 6 . 2 ) 


We now show that for N an odd integer, the reciprocals of the roots of the 
polynomial 


C(D) = Cq + c^D + + ••• + c^d'^ (6.3) 

are the desired mass points, x^, X 2 , **', and the probability masses at these 

points, to^, ^ 2 , •••, given by (3.15), do indeed yield the approximate 
probability 


Pr(X = x^) = 0 )^ ; £ = 1,2, •••, v (6.4) 

which is the unique solution to the moment problem given moments p^, ••*, p^^. 

In the following, if X is a discrete random variable, we assume that the 
true probability distribution has at least v points with nonzero probability. 
Otherwise there would be no point in finding an approximating probability 


*Except for pathological cases, we have v = 


N + 1 

2 • 
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distribution for X, Now note that since = 1, (6.1) can be expressed in 
matrix form as follows: 



\)+l 


2v-l 


'v-l 


= 0 


(6.5) 


This corresponds to v linear equations in v variables c^, c^, and has 

a unique real solution if 


^0 

^1 ^2 


v-1 

% 


M = 


( 6 . 6 ) 




'2v-2 


is nonsingular. M is singular if and only if there exists a column vector ^ 
with elements a^, a^, ••*5 that 


a"^Ma = 0 


(6.7) 
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or 


v-1 v-1 



i=0 j=0 


a.a.y , . . 
1 J i+3 


= 0 


Recalling the definition of the moments, then equivalently 


( 6 . 8 ) 


i-0 j-0 


or 


E 




- 0 


( 6 . 10 ) 


This is possible only if all the values of t ic random variable X are at the 
v-1 roots of the polynomial 


A(x) 


v-1 

i 

a.x 
1 

i=0 



( 6 . 11 ) 


For our case, this is not true since we assumed that at least v points have 
nonzero probabilities. Hence, M is nonsingular and the Berlekamp-Massey 
algorithm yields a unique solution given by the polynomial C(D) in (6.2). 

The roots of the polynomial C(D) must be distinct and real. To see this, 
we consider the reciprocal polynomial 

Q(D) = 

2 V 

= c + c -D + c + ••• + c^D (6.12) 

V v~l v-2 0 
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and show that the roots of Q(D) are distinct and real. 

Suppose, for m < v, X 2 , are the only distinct roots of Q(D). 

Let ^ 2 * *•*» (ro' ^ m) be those real distinct roots where Q(D) changes 

sign for real D. Define polynomial 

m' 

R(D) = ]^ (D - 
i=l 

( 6 . 13 ) 

i=0 


Then 


Q(D)R(D) ^ 0 (6.14) 

for all real D since changes in sign of Q(D) are reversed by sign changes in 
R(D) . Also the only real numbers for which 

Q(D)R(D) = 0 (6.15) 

are the roots A-, •••, X . Since the random variable X takes on values at 

X / m 

other points besides these m root points (m < v) we have 

E[q(X)R(X)] > 0 (6.16) 
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However 


E[Q(X)R(X)] = E 


m' 




L\j=o 


^i=0 


= E 


m / V 

E^(E^. 

|_i=0 \j=0 


X 


V+i-j 


m' /V 


E^.Ev..- 


i=0 \j=0 


which equals zero since 


V 

E c.p . = 0 for n > 

J n-3 - 


j=0 


Thus, by contradiction, we must have m = v and all roots of Q(D) and C(D) 
real and distinct. 


Since all roots are real and distinct 








must be real since the polynomials P(D) and C(D) have real coefficients. 

also know that co- , , •••, to satisfies 

12 V 


V 

E 


’ k = 0.1.2, N 

£=1 
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tJ 


(6.17) 


(6.18) 


must be 


(6.19) 


We 


(6.20) 


1 



Hence, for any polynomial F(X) of degree ^ N, we have 


V 

e[f(X)] = ^ 




F(x,) 




( 6 . 21 ) 


Choose, for some 1 < _< v. 


F(X) = Y\ ° 

j=l 


which has degree 2v - 2 _< N. Then 


V 

E[F(X)] 

i=l 


( 6 . 22 ) 


V 

j=i 


> 0 (6.23) 


and thus 


0) > 0 (6.24) 

Xj 


The condition 



(6.25) 
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completes the proof that is a set of discrete probability 

weights. 

It is also easy to see that this set is unique from the constraints of the 
moments given by (6.20). In matrix form, this is 


1 

i 


1 

1 

1 .... 

— 

1 


1 

3 

1 



^1 

^2 

X3 .... 

X 

V 


“2 

^2 


2 

^1 

2 

^2 

2 

X3 .... 

2 

^v 


• 

• 


• 



• 


• 

• 


• 



4 


• 

. 

I 

• 



« 


• 

• 


• 



• 


: • 

^-1 


v-1 

v-1 

Xz 

v-1 

• • • 1 

V -1 

X 

V 


‘"v 

■— — 







^ — 


(6.26) 


where the matrix 


M = 


v-1 v-1 v-1 

Xf x^ x^ 


V 

2 


v-1 


(6.27) 


is a Vandermonde matrix (Ref. 11) with nonzero determinant since x^^, 
are distinct. 
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VII. Generalization to Correlated Random Variables 


In many communication systems such as the satellite transponder system 
example in Section I, we want to evaluate the expected value of a function of a 
complex random variable such as 


W = X + jY ; j = /T (7.1) 

This is typically the complex envelope of a narrowband signal. If we follow our 
earlier approach and assume we have available a set of complex moments 

= E(wS ; k = 0,1,2, •••, N (7.2) 

then we can again apply the Berlekamp-Massey algorithm. The Berlekamp-Massey 
algorithm works for any field so certainly the complex number field is no 
problem. This algorithm, in fact, was originally developed for finite fields. 

Despite what seems like an obvious extension of the previous results, the 
complex random variable generalization of the moment technique needs to be 
investigated further as there are some special cases where it does not seem to 
work. Suppose, for example, 


W = 


(7.3) 


where 0 is uniformly distributed over (0,2 ir). Here, we have 


which yields the trivial uninteresting solution 


Pr(W = w) = 


II 

0 

; k = 

1,2,3 

solution 


r ‘ 

o 

II 

(o ; 

w 0 




(7.4) 


(7.5) 
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In this case, however, we can reformulate the basic desired expectation as 


E[G(W)] = E[G(e^®)] 

= E[H(0)] (7.6) 

where we achieve an approximation using moments of the real random variable 0. 

Another example which causes problems is W = X + jY where X and Y are 

2 

independent zero mean Gaussian random variables with variance a . Then 

W = (7.7) 


where A is a Rayleigh random variable that is independent of 0 , a uniformly 
distributed phase random variable. Again we have complex moments given by (7.4) 
yielding the trivial approximation (7,5). This case can also be solved easily 
using a reformulation as follows: 

E[G(W)] = e[G(X + jY)] 

= e[f(X,Y)] (7.8) 

Now we can apply the real random variable approximation for X and Y to 

obtain 


fr(X = = Pr(Y = y^) = = 1,2, •••, v 


based on moments E(x'^) = E(Y^); k = 0,1,2, ••*, N. Then 


(7.9) 


E[F(X,Y)] 


EE 

£=1 m=l 




m 


(7.10) 


This, in fact, is the double application of the Gauss-Quadrature rule for 
Gaussian integrals called Gauss-Hermite approximation. 
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Tht^ above pathological cases can be easily handled using the single real 
random variable moment technique described in Sections II through VI. They 
point out , however , the need to further investigate the complex random variable 
generalization . 

We now consider the generalization to two correlated real random variables 
which includes the complex variable problem as a special case. As shown next, 
this approach requires multiple application of the single random variable tech- 
nique and, most Importantly, does not result in unique solutions . 

Assume we wish to evaluate E[F(X,Y)] when we only know the (N + 1) joint 
moments 


= E(xV) ; i, k = 0,1,2, •••, N (7.11) 

2 

We issume the joint probability of X and Y is approximated by v pairs of 
po-'nts* 


^VmU^ ’ m = 1,2, V 

and probability masses at these points given by 

fr(X . xj, Y . ; H, m - 1,2 , • ■ • . x (7.12) 

where is the approximation to the conditional probability of Y = Y^|£ given 

X = while 0 )^ = Pr(X = • 

This allows for the approximation 


e[f(x,y)] = 


V V 
£=1 m=l 


n, I 


(7.13) 


*The notation indicates that the discrete set of points at which Y will be 

allowed to have probability mass depends on the discrete set of points chosen 
for X to be allowed to have probability mass. 
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To find the approximating joint discrete distribution, consider first the 
constraints imposed by the given joint moments of (7.11), namely. 


= E(xS^ 

V V 

* 1, k - 0.1,2, •••, B 

£=1 m=l 

For k = 0, we have 

V 

>'10 ^ “ 

£=1 


(7.14) 


(7.15) 


By applying the single real random variable moment technique, we can find the 
smallest set of v [v = (N + l)/2 for N odd except for pathological cases] unique 
mass poiiyfs x^^, and weights ojj^, ( 1 ) 2 , • • • , (o^ satisfying the N + 1 

moments of (7.15). 

Next observe that if we define the approximating conditional moments* 

V 

“ XI ’ k = 0,1,2, ..., N (7.16) 

m=l 


*Note that this approach does not insure that the approximating conditional 
moments be equal to the true conditional moments Pic| = E(Y^Ix = x^) nor 

does it guarantee that they are a valid set of moments in the sense of producing 
a convergent Berlekamp-Massey algorithm. More often than not, however, the 
approach will be successful and yield meaningful results. 
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for each Z = 1,2, •»*, v, then we merely apply the real random variable moment 
technique v times to find the points conditional probabilities 

{Pmj^)' We get these conditional moments from the following “'fpresslon; 




For each fixed Ic we have a set of linear equations for 

since {to^} and are known. These equations can be expressed in the matrix 

form 




(0 

V 



(7.18) 
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where the transformation matrix is M of (6.27). Defining the LaGrange 
polynomials (Ref. 11) 



n = 1,2, 


V 


= ag(n) + a^(n)D + 02 (n)D^ + ••• + a^_^(n)D^ ^ 


with the obvious property that 


T 

n 


(Kj) 



^ = n 

JZ. ^ n 


(7.19) 


( 7 . 20 ) 


then, equivalently the coefficients ot.(n); n = 0,1,2, •••, v - 1, which are 

* ^ 
easily found, have the inner product property 


jaQ(n), a^(n), •••, a^_^(n)J 


2 


X. 


v-1 



il = n 

(7.21) 

Z ^ n 


*See Appendix A. 
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Hence, the inverse of M in (6.27) is easily seen to be 


M 


-1 


a^d) ot^(l) 062 ( 1 ) 
Oq( 2 ) a^( 2 ) 02 ( 2 ) 


oij^(v) a2(v) 


Vi‘^> 


a , (v) 
v-1 


(7.22) 


and from (7.18), the desired conditional moments are found from the joint moments 
by 


“l \ll 


>^0k 

y\ 



CM 

CM 

3 


^Ik 

• 


• 

• 

= M 

• 

03 y, 1 
V k V 


Vi,k 



— 


; k = 0,1,2, N 


(7.23) 


This completes the solution. ’ 
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The above solution to the two correlated real random variables moment 
problem is clearly not unique since we could have interchanged the role of the 
random variables X and Y. Also it is not clear if this approach results in the 
fewest number of mass points compatible with the given joint moments. Finally 
this procedure may be improved by using an invertible transformation ^ to define 
new variables X and Y where 



(7.24) 


^ i Ic 

Joint moments = E(X Y ) can be easily found from the original moments and we 
can easily find F(*,*) such that 


e[f(x,y)] = e[f(x,y)] 


(7.25) 


The choice of new transformed variables would come from examination of the original 

physical problem that led to the requirement for evaluating E[f(X,Y)]. Joint 

moments may al^o be easier to find by an appropriate transformation. For special 

cases such as when X and Y are correlated Gaussian random variables, we can 

always find a transformation such that X and Y are independent zero mean Gaussian 

2 

random variables with variance a =1. Then the problem of evaluating the expec- 
tation of F(X,Y) = F(X,Y) reduces to a double application of the single real 
variable moment solution. 

An alternate approach to the correlated random variable problem, which 
does not depend on whether X or Y is chosen as the unconditioned random variable, 
is based on a direct two-dimensional generalization of the one-dimensional solu- 
tion. In particular, we do not search for pairs of points and associated proba- 
bility masses whose values for the second dimension are conditioned on those 
found for the first dimension. Rather, we directly proceed to find joint mass 
points 

; l = 1,2, ‘ V (7.26) 
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and weights u ; £. = 1,2 v at these points giving the approximate joint 

probability distribution 

Pr(X = x^, Y = y^) = ; £ = 1,2, •••, v (7.27) 

2 

Here again the available input consists only of the (N + 1) joint moments of 
(7.11) and the desired output is the evaluation of E[F(X,Y)]. Once we have the 
approximate joint probability of (7.27) we may make the approximation 

V 

E[F(X,Y)] s; y^^g),F(x^,y,) 

A=1 


Our goal is to find an approximate joint probability distribution as given 
in (7.27) with the fewest nmber of points v that satisfy the joint moment 
condition 


= E(xV^) 


- L 


i k 


i ,k — 0,1, 


N 


(7.28) 


£=1 


First, we denote x , x , • • • , x^ as 

X ^ X 

set x^, x^, •••, x^. Similarly, we let y^, 
bers among the set y^, y^, •••, y^^. Thus, 
pairs (x^, y^) and the desired set of mass 
subset of all such distinct pairs. 

Next, define the polynomial 


the set of distinct numbers among the 



there are a total of v v < v distinct 

X y 

points (x^, y^); i = 1,2 •••, vis a 


<=x®> ■ 


A 

n 

Z=1 


(1 - 




= Sq + a^D + • • • + D 

X 


(7.29) 
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where = 1. Note that analogous to (3.6), 


X /V 


W'V 


m=0 


m =0 \ 1=1 


X 


= L 


a X. 
m i 


-Til 


il=l 


m=0 


^ V<=x ( 


£=i 


= 0 


(7.30) 


-1 .V -1 

since x„ (or x ) is a rooc of C__(D) 
m A 

alternate form 


Thus, (7.30) can be written in the 



(7.31) 


For a given value of j, (7.17) has a shift register interpretation analogous to 
Figure 3. In particular, the conditions on the coefficients a^, 82 a-^^ 

imposed by (7.17) are those of a v -tap feedback register that is required to 
be able to generate N + 1 different sequences, namely. 
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^ 10 > % -1,0‘ ^ .0’ ■'^NO 

X X 

"ll’ •••• "v -1,1’ ^ 

X X 

-1,N’ ,N’ ^NN 

y— 

Initial Condition 

For each of the above N + 1 sequences of length N + 1, the first v terms serve 

X 

as the initial loading of the feedback shift register specified by the polynomial 
coefficients a^, •••, a^ • 

X 

Next define 

V 

y 

£=1 

V 

= bp + b^Z + • •• + Z y (7.32) 

where b^ = 1. Using the same development as that leading to (7.30), we obtain 
now 


^ 00 ’ 

^01’ 

'^ON’ 
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V 

y 


E 

n=0 


by.. 


V 

y 



n=0 






^ / ^y 



o)„x„Vo^c^(y 


“r£ 


z=i 


z 


= 0 


( 7 . 33 ) 


since y,-l (or y is a root of C„(Z). Thus, analagous to (7.31), we can write 
^ m I 

(7.33) in the alternate form 



(7.34) 


Note that, although at this point, we do not the two sets of points 

x^, •*"> Scvx 72 » ***^ yvy their gene v-'- I polynomials 

Cy(Z), we have the interpretation that the given (N ^ V) int moments of 
(7.11) are generated by two linear feedback drift regxcnx' with feedback tap 
coefficients that specify these polynomials. This is a new interpretation or 
formulation of the classical two-dimensional moment problem. Furthermore, even 
after we were to find these two sets of points by some suitable algorithm, they 
would not yet be paired together. Thus, at that point, it would still be unclear 
which V pairs (x^ , y.); - lj2, •••, v out of the v v pairs (x. , y ) are 

valid mass points. 
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To resolve this ambiguity, we proceed as in the one-dimensional case by 
next defining the polynomial 


V X 


y r 


p(D,z) = ~ 


Z=1 j=l 


1=1 


V -1 V -1 
X y 



I / 


(7.35) 


i=0 ; j=0 


where the primes on the two products in (7.35) respectively denote omission of 
the factors corresponding to x^ = x^ and 9^ ~ 7 Thus, we may write the equiva- 


lent relations 


(1 - ^Xj) 


n 


n'a-^V- 


i=l 


n<i - "V 


.1^ 


(1 - zy,) 


(7.36) 


Also define the joint moment generating polynomial 


y(D,z) 


EE- 

i=0 j=0 


IJ 


.D^Z-^ 


Then, assuming the moment relationship of (7.28), we get 


(7.37) 
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J 


00 00 / \) 


1=0 j=0 \ A=1 




V 

E' 

ji=i 


E 

1=0 


<v 


L j=o 


V 

= E‘ 

a=i 


'Z (1 - DXj^) (1 - Zy^) 


(7.38) 


which when multiplied by Cj^(D) and C^(Z) produces the relation 


V 

y(D,Z)C^(D)CY(Z) = ^ 


to, 


£=1 


X 


n 


(1 - Ds;) 


(1 - DXjj^) 


V 

y 

1 =;. 


n ■ "^ 1 * 


(1 - zy,> 


= P(D,Z) 


(7.39) 


Equating coefficients of equal powers of D and Z In (7.39) yields the coefficients 
p. , of the polynomial P(D,Z). The procedure for accomplishing this is as follows. 


Substitute the polynomial representations of y(D,Z)» C^(D), and C^(Z) 

A X 

given In (7.37), (7.29), and (7.32) respectively into the product in (7.39) to 
yield 


p(D,Z)Cjj(D)Cy(Z) = 


V V 

X y 00 00 

EEEEW*'” 

k=0 Ji=0 1=0 j=0 


V V 

X y 00 00 



k=0 £=0 i=k j=£ 


(7.40) 
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Note that 



min(j ,Vy) 


< >^E E ,< > 

£=0 j=£ j=0 £~0 


X oo 


min(i,v ) 


' >-E E < > 

k=0 i=k i=0 k=0 


Using the equivalences of (7.41) in (7.40) yields 


00 00 


min ( i > v^) min ( j , ) 


y(D,Z)Cjj(D)CY(Z) 




i=0 j=0 k=0 £=0 


V -Iv -1 
X y 


3-.7J 


= P(D,2) = 2^ 2_j 

i=0 j =0 


Finally , 


1 2 

'’M ■ EE “kVi-k.j-<! • i 1 

i =0 £=0 

j = 0 , 1 , •••, - 1 


This relation is Che two “-dimensional generalization of (3.11) 


Note that for i > v^, we have from (7.42) the condition 


X 


]^®k'^i-k,j-£ ° 


k=0 


(7.41) 


(7.42) 


(7.43) 


(7.44) 
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and for j > V , the condition 

J y’ 



^)l^i-k, j-S, 


= 0 


(7.45) 


both of which agree with the conditions on a^^, and b^, •••» 

previously found in (7.31) and (7.34) respectively. 

In summary, given the polynomials Cjj(D), Cy(Z) and the known joint moments 
of X and Y, we can easily obtain the polynomial P(D,Z). Given C^(D), 0^(2), and 
P(D,Z), we show next how the weights u\^; £ = 1,2, •••, v are found. 

Prom the definition of C (D) in (7.29), we have 

A 

" 35 V’» 



(7.46) 


— I 

Thus, to evaluate (x^ ) where x^^ = x„,q, only the term corresponding to m = m^ 

in the summation of (7.46) would have a nonzero contribution, i„e.. 


V 
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where the prime is again used to denote omission of the factor in the product 
for which x. = x . Similarly, for C (Z) , we would have 

J JC A 


k") - 


® V”> 


D = y 


-1 


m 


= - y 


m 


j 

n 

i=l 


'/ y.,- 


'my 


(7.48) 


From the definition of P(D,Z) in (7.35), we observe that 


(7.49) 


Also , 


- rt'(' 4) a 


1 - 




or 





\ 

K ’^1 ) . , 


T 5 X/ — ,V 

' ^y( yji" ) 


(7.50) 


(7.51) 


The above relation for the weights of the approximating joint probability dis- 
tribution is clearly seen to be the two-dimensional generalization of (3.15). 
Also, we have demonstrated that out of the total of pairs of points 

(x. , y ) only V of these pairs, namely (x , y. ) ; ^ = 1>2, V will result in 

nonzero probability weights as determined from (7.51). 
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As a check on our procedure, let us examine the known special case where 
X and Y are independent. Here the joint moments have the form 

^ij ^ •••, N (7.52) 


where 


= E(x"-) ; 1 = 0,1, •••, N 
(yY)j = E(Y^) ; .i = 0,1, •••, N 


(7.53) 


Here (7.31) reduces to 


(Viy) , = 

A 1 


.A. 

-E 

m=l 


m A 1 . 


•m 


X ^ V 


X 


(7.54) 


which is the single sequence shift register requirement as in (3.7). Similarly 
(7 . 34) reduces to 


(Py^j 


E 

n=l 


''n'Vj-n 


3 s v 


(7.55) 


Thus, we get the correct sets of points x^^, x^j and y^, ^ 2 ^ •••j Yy 

where i.e., all pairs (x^, y^) have nonzero probability weigtits. 

Suppose now that (wx)mj " 1»2, •••, and (w^)^; ti == 1,2, •••, v are the 
probability weights for each random variable. Then, (7.38) has the form 
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y(D,z) = 


(7.56) 


V 

X 


V 

y 




("'v)™("’v) 




m=l n=l 


1 

(1 - DS_„) ( 1 - zyj 


and corresp j)ndingly (7.35) becomes 



(7.57) 


(7.58) 


Thus, in conclusion, we see that the general two-dimensional forms of the results 
are consistent with the known case where X and Y are independent. 


VIII. Constrained Moment Problem 

In some applications, we may wish to place a constraint on the mass points 
when solving the moment problem. In this section, we consider a few special 
cases where some of the probability mass points are fixed. Here let X be a 
random variable with moments 


- E{X"^} ; k = 0,1,2, - , N 


We want to find an approximate discrete distribution for this random variable 
based only on the given moments. Suppose, however, we require that the approxi- 
mate distribution have probability mass at given points y^, •••, y^. Our 

goal is to find the fewest points x^, x^, •••, and probabilities 


Pr(X-xp=o)^ ; i=l,2, 


Pr(X = yj) = 2 ;^ ; j = 1,2, • • • , p 


that yields 




k = 0,1,2, 


Hence, given moments y^, y^, •••, y^ and a set of fixed points y^, y^, ‘‘’j y^j 
we wish to find the smallest v, mass points x^, x^, ’•*> and probabilities 


^ 1 ^ ^ 2 ^ ***> ^ ^ 1 ’ ^ 2 ^ 


, z where 
P 


E-.-S 


z =1 
p 


r/r-E ei-A*--" j:ct 
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Suppose the unconstrained moment problem yielded V* and x*, x*, • • • > 

0 )*, '*)*, - ' • , Let a,b be any numbers where 

a < min x* 

" 1 ^ 

(8.5) 

b > max X* 

- i ^ 

We noxiT examine some special cases of the constrained moment problem [see Krein 
(Ref. 1), pp. 53-55]. 

Case I: p = 1, y^^ = a 


i=l 


k = 0,1,2, • • • , N 


( 8 . 6 ) 


Consider 


- ay 


k+1 “^k \i=l 


> a>iXi +z^a ! - a^ o).x. + 
\i=l / \i=l 


p, - a 


V 

E W.X. + Z„ 
11 1 

^ i=l 


a I ~ a 




a) 


i-1 


V 


12 "I'i - 


(1 - z^)a 


i=l 

V 


12 Vi'*i - 


1^1 


E-. 


io . (x. - a) 


j=l 


(8.7) 
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Defining 


V 

since 1 - 

i=l 


0 ) . (x . ~ a) 
11 ' 


1 V 


i = 1,2, •••, V 




a) 


j=l 


( 8 . 8 ) 


we see that 


and 


(3. > 0 since x. > a (8.9) 

1 — 1 — 


V 



i=l 


( 8 . 10 ) 


(8.9) and (8.10) reveal that the set of weights i ~ 1>2, •••, v} has the 

properties of a probability distribution. Substituting (8.8) into (8.7) gives 


y 


k+1 


a 



k = 0,1,2, ■ • • , N-1 


( 8 . 11 ) 


Hence, given •••, y^, compute new moments. 


k+1 


a 


k = 0,1,2, • • • , N-1 


( 8 . 12 ) 
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and use the unconstrained moment solution to get X 2 » and 

Note that 




1 



- a) 


so that 


A 

OJ . 
1 


0 ) . (x. - a) 

X 1 
Pi “ a 


( 8 . 13 ) 


J 


( 8 . 14 ) 


or 


m.(Pl - a) 

03 , 

1 X. - a 

X 


i = 1,2, 


V 


( 8 . 15 ) 


and 



( 8 . 16 ) 


We thus find the solution to the moment problem where -(8. 6) is satisfied with 
one mass point = a fixed and all other mass points having values greater than 
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Case II: p = 1, = b 


Consider 




V 

= (o^x^ + ; k = 0,1,2, • • • , N 


i=l 


(8.17) 


- ^k+1 
b - y. 



\ 0) , (b - X . ) 


1 1 


i=l 


V 

E 

j=i 


0) . (b - X. ) 
2 2 


(8.18) 


V 

since 1 - . Defining 


j=l 


u . = 


w . (b - X.) 
1 1 


1 V 


; i = 1,2, ••• V 


E w.(b - X.; 

J 1 


j=l 


(8.19) 


we see that 


CO. > 0 since b > x. 
1 — — 1 


( 8 . 20 ) 
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Hence, given y^, •••, y^^, compute new moments 



k = 0,1,2, • • • , N~1 


and use the unconstrained moment solution to get x^, x^, • * • , x^ and 
Note that 



so that 






Case III: p = 2, y^ = a, y 2 = b 




I 

4- z^a + z^b 


i=l 


Consider 



V 

^ (”^i (a + b)x^ ~ abj 

i=l 


4- z 






- X )(x - a) 




Thus, 


-“k+z + '>>''k+l - “'’''k 1-1 


V 

- ’tiHki - k) 


-^2 + (a + b)vij^ - ab 


E CO (b - X ) (x - a) 
3 3 3 


j=l 


(8.30) 


Define the new probability distribution 


10 . = 


10, (b - X.) (x. - a) 
X X X 


i V 


i = 1,2, • • • , V 




(8.31) 


and moments 


— -\-k2 + '■>'^k-H - ^’’“k 

*^k -p„ + (ct + fa)u, - ab 


k = 0,1, • . • , N-2 (8.32) 


Thus, using moments Pq, •••, find x^^, Xg, •*•, x^ and 

tio^, toj, • • • , (0^ from the unconstrained solution. Then 


w . 

1 

-P 2 + (a + b)p^ - ab 

** • T T 0 a M A ^ \ 

(jl) . — 
1 

(b - x^) (x^ - a) 

9 X — X,iC, 


(8.33) 


To find and Z 2 solve 


V 

= X) “i ^ 


i=l 

V 


It = 7 + Z.a ~\r 2„b 

i jL^ X X 1 2 

i=l 


(8.34) 
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IX. Accuracy o f the Moment Approximation 


In this section, we examine the accuracy of the moment approximation for 

Ic 

E{f(X)} where X is the random variable with given moments ” E(X ); 
k = 0,1, ***, N, The solution to the moment problem yields points {x.} and 

weights where we have the approximation 


Pr(X « x^) = ^ 


(9.1) 


and 


V 

Jl=l 


(9.2) 


Tx\to types of bounds are oresented for the accuracy of this approximation. The 
first bound assumes a bounded K + 1 st^ derivative of f(x) while the second bound 
assumes that X is a bounded random variable and the N + 1 _st derivative of f(x) 
is convex 0 or convex U in the finite range of X. 

A. Bounded Derivative 

Assume all K + 1 derivatives of f(x) exist everywhere and that 


j(K+l) 4 ^ 


K+1 


dx 


K+1 


f(x) 


(9.3) 


is bounded for all x. That is 

If^^'^^hx)! for all x 


(9.4) 


*For N even we take K = N - 1 while for N odd we take K = N - 2. 
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Next, consider integration by parts to obtain 


'X 




u) 


n-i 


du 


(n ~ 1)! (n - 1)! 


■t 




JO 




n-1 


(ti - 2) ! 


du 


n-1 


k 


+ f(x) 


k=0 


Thus, setting n = K + 1, this becomes 


(9.5) 


K y, r-> 

f(x) =Y^ fr+ I f'‘^'"^(u) du 


:(K+1). S (x - u) 
k! ■ I " K! 


K 


k=0 


Jo 


(9.6) 


and changing the variable of integration, 


K 


.K 


f(x) = f (0) "j- + I f(K+l) 


f^‘"""^xu)du 


k=0 


'0 


(9.7) 
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Now using the bound on in (9.4), we have 


1 1 I = 

' K+1 ' 


^K+l r (1 . u)*^ 

^ 1 K'M| 
1 |x 1 

r (1 -_u)" ((K+1) 

Jo 


Jo 


0 


.. -D I K+1 I I (1 - U)^ , 

- \+l I I ’ K! 

'o 


K+1 1 


B 


X 


K+1 (K + 1) ! 


(9.8) 


Since the first N moments are the same for the true and approximate 
probability distributions, we have 


/ k=0 


K 


= E 


\jLj 

k=0 




(9.9) 


Thus, the approximation error is due to tae integral cerm in (9.7). Taking the 
expected value of (9.7) with respect to the true and approximating distributions 
and differencing the results yields the error bound 

^ |E{f(X)} - E{f(X)}| 

I^Uk+i} - £ |e{Ij^_^^}| + |e{Ij^_j_j^} I 


< B 


K+1 


e{|x'^«|} + e{|x"«|} 

(K + 1) ! 


(9.10) 
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J 


If, as assumed, K = - 1 for N even and K - N ~ 2 for N odd, then 

~ £qj- even or odd 

and 

B{i/«|} = e{|/«|} . 

Thus, substituting (9.12) into (9.10) gives the desired result 


(9.13) 

The bound derived above can be generalized to functions of two correlated random 
variables. 

B. Bounded Random Variables 

Suppose X is a bounded random variable where 

a£X£b (9.14) 

Given moments • • • , p^, we now define principal probability distribution 

functions. These are approximate probability distribution functions subject to 
various constraints on mass location points of the type previously considered in 
Section VIII. 

Case I: N = 2n - 1 (N odd) 

For this case, the principal* distribution functions are the solutions to 
the following: 

(1) Unconstrained; V = n 


Xi, 

x^. 


0)^, 

“ 2 ' ’ 

* * , 03 


96 




(9.11) 


(9.12) 


1 



(2) Constrained: p = 2, Yt == a, y« = b, v = n 1 


* ^2 ^ 


> ^ 1 
n-1 


• • • , W. 


n-1 


and 


52l, Z2 


Case II: N = 2n (N even) 

For this case, the principal distribution functions are solutions to the 
following: 

(1) Constrained: p- 1, y^-a, v~n 

x^, 


0)- , 03^ , • • * , a). 


n 


and 


(2) Constrained: p“l, y^==b, v = n 


^1 ^ ^2 ^ 




• • • 5 


and 


Denote the two principal distribution functions as 


S «- *u> - “u - * ■ "i 


(9,15) 


and 


Pr2 (X = =0)2^ ; 1 = l,2y ’ ” , ^2 


(9.16) 
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An important bound due to Krein (Ref. 1) is given as follows: 
Let f(x) be any function where 


f(N+l) 


(x) 


,N+1 


dx 


N+1 


f(x) 


(9.17) 


is either convex U or convex 0 in [a,b]. Then 


Ej^{f(X)} < E{f(X)} < E2{f(X)} 


(9.18) 


where 


E^{f(X)} = 


£=1 


(9.19) 


and 





(9.20) 


Note that if 


f(N+3)(^) Q ^ ^ 


(9.21) 


then f 


(N+1) 


(x) is convex U In [a,b] whereas if 


f(N+3) 


(x) < 0 


for all X E [a,b] 


(9.22) 


98 



then is convex f) in [a,b]. Yao and Biglieri (Ref» 3) have applied 

these results to the Gaussian probability integral f(x) = Q(x) [see (2.15)] 
to obtain tight bounds on error probability performance of BPSK signaling over 
additive ^hite Gaussian noise channels with bounded interference signals. 


99 



X. 


Conclusions and Other Applications 


Our primary motivation for this study of computational techniques based on 
moments is the evaluation of satellite communication system performance with 
uplink interference signals and satellite nonlinearities. Here we presented new 
ways of solving the moment problem, examined the accuracies of the approximations, 
and extended the techniques to two correlated random variables. The computational 
requirements are modest and the approximations are very accurate for evaluating 
bit error probabilities (Ref. 3). 

Although our example stressed evaluation of bit error probabilities, we 
can apply these moment techniques to the evaluation of other parameters, such as 
channel coding cutoff rates under both normal and mismatched receiver cases 
(Ref. 12). Most modulations and interference signals can be handled using these 
moment techniques . 

Another very important application of the computational techniques based 
on moments is the determination of the probability distributions of the outputs 
of a discrete-time dynamical system. Specifically, consider a discrete-time 
system with inputs that are independent random variables with known probability 
distributions. Figure 6 shows a generic system where has a known probability 
distribution and are independent random variables with known probability 

distributions. There are many examples of control systems, queueing systems, 
and synchronization systems where this type of model occurs. Our goal is to 
find approximate probability distributions for the state X^^ at time, t^, 
k = 1,2, i.e., we wish to determine an approximate probability distribution 

of the form 


■ “k* 


£ == 1,2, • •• , V 


( 10 . 1 ) 


for k = 1,2 , • * 



Figure 6. Discrete-Time System 
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The moment technique can be used in a recursive manner to solve this 
problem as follows: 


Step 1: 


Compute the moments of 


''ik’ ® 


w 

- ”o>) * 


k = 0,1,2, N 


( 10 . 2 ) 


where we use E(*) to denote expectation over both the initial condition random 
variable Xq and the input variable n^. 

Step 2: Solve the moment problem to obtain the approximation 


Pr(X^ = = 0 )^^ ; a = 1,2, •••, v 


(10.3) 


Step 3: Compute the approximate moments of X 2 using the probability distribu- 

tion obtained in Step 2 for computing the expectation over X^, viz., 



V 


(10.4) 


where E(») now denotes only the expectation over the variable n^. 

Step 4: Solve the moment problem to obtain the next approximation 

Pr(X 2 = x^^) = ; )!• = 1,2, •••, v (10.5) 
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F. 


etc. By repeating this procedure, we obtain 

" \l ’ S, =• 1 , 2 , • • • , V 

k = 1,2, • • • , (10.6) 

as desired. Since in each step we use valid moments, the algorithms for solving 
the moment problem should ’*ot encounter any difficulties. Increased accuracy 
can be achieved by increase g lha number of moments used in each stage. Indeed, 
we could consider using different values of N at each stage. 

Note that the above procedure does not require that the system, which is a 
Markov process, be irreducible. Also, we can extend the results to second order 
processes of the form 


■ k-0,1. ••• 


( 10 . 7 ) 


Here we can define 


Y 


k 


= X 


k-1 


and obtain the vector first order form 


( 10 . 8 ) 


\+l 



Y 

— 


__ k+1 


L K J 


( 10 . 9 ) 


{ 

This is a special case of two dimensional systems of the form 


\+i ~ 

( 10 . 10 ) 

\+l “ 

where Xq>^q tiave known joint probability distributions and ^ 

sequence of independent pairs of random variables with known joint probability 
distributions . 
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To find an approximate joint distribution for (X, ,Y ) of the form 

■ *kr h • ^kmlll) ‘ Pkm|n“kt ■ (W-IO 

for k = 1,2, •••we can repeat the steps given above using the joint moments and 
the two random variable generalization of the moment technique discussed in 
Section VII. 

Here we have demonstrated an application of the computational techniques 
ba^ed on moments to two dimensional first order Markov processes. Many special 
cases of this application need to be further explored. Synchronization systems, 
in particular digital phase-locked loops, fall nicely into this category. 
Queueing systems analysis is another area where such techniques ^^’ill be very 
useful . 

The computational evaluation technique based on moments presented in this 
report is a very general and powerful numerical technique for evaluating the 
performance of a wide range of systems particularly communication systems. We 
feel that the applications of these moment techniques have just begun. Subse- 
quent reports will be devoted to the analysis of various modulation and coding 
schemes used over satellite channels where the techniques described here will be 
the basic analytical tool- 
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APPENDIX A 


A Recursive Method for Finding the Coefficients of a 
Polynomial Generated by a Product of 
First Degree Factors 


Consider first the problem of determine the coefficients {a^(v)} of the 
polynomial 


P„(D) -I~[ O’ - V 

«.= 1 

= a^Cv) + a^(v)D 4- a 2 (v)D^ 4- • • • 4- a^(v)o'' (A-1) 


We start by defining 


P^(D) = D - = a^d) 4- a^(l)D 


(A-2) 


Thus , 


aQd) 


-X 


1 


a^d) = 1 


(A-3) 


Next, consider 

P 2 (D) = (D - Xj^) (D - X 2 ) = P^(D)(D - X 2 ) = SqCZ) 4- a^(2)D 4- a2(2)D^ (A~4) 

Clearly then 


aQ(2) = -X2(-x^) = -X2aQ(l) 

aj^(2) = -X 2 (l) + (l)(-x^) = -X 2 a^(l) 4- agCl) 

a2<2) = d)(l) = a^(l) (A-5) 


"F.DING PA'GE BLAKIC NOT FILMED 
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Generalizing to arbitrary k, we define 


k+1 

Uj.1 

^ 4- 1) + (k + 1)D 4“ ^2^^ ^ -f • • • ^k4-l (A-6) 


and hence 


a^Ck + 1) = -Xj^^^apCk) 
3j^(k + 1) = ■*■ 


a£(k + 1) " ^ a^(k) 


a^(k + 1) = -x^+iaj^(k) + aj^_^(k) 


a^^ldc + 1) = a^(k) 


(A-7) 


Finally, letting k = v ~ 1 in (A-7) gives the desired result, namely a recursive 
relation for the coefficients of the polynomial in (A-1) , 




Now referring to ^7.19), we are interested in determining the coefficients 


{a^ '^(v)} of the polynomial 


Q®)(D) ® - V ■ “ 

a=i 


P,.(D) 


D - X 


n 


a^'^^(v) + aJ"^(v)D + a^’^^(v)D^ + * 


+ a ' (v)D 
v-1^ 


n = 1 ,2 , • • • , V 
(A-8) 
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The procedure to be followed is Identical to that used in the root-finding 
algorithm associated with the Berlekamp-Massey algorithm discussed in 
Section IV. There a recursive procedure was described for removing a first 
degree fa .:or from a known polynomial to arrive at the coefficients of the 
reduced polynomial. Applying that procedure to this case results in 

aj^(v) = -x^aj^^(v) + ay’^^(v) 

a2(v) = 


or equivalently, 




1^. 


a T (v) = -X a^^T (v) + a^'^hv) 
v-1' n v-1 ^ v-2^ ^ 


a (v) = (v) 

v-1 


(A-9) 


a'^'cv) - 


aj^(v) - Sq^'cv) 






a'"?(v) = 

V-1 


- a<”>(v) 


X 


-• a (v) 

V 


n = 1,2, • • • , V 


(A-10) 
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A special case occurs if = 0 for any n. In that situation (A-10) is replaced 
by 


(v) = a^ (v) 
a{“^v) = a^(v) 



(V) 


av(v) 


(A-11) 


Finally, comparing (D) with T^(D) of (7-19), we immed'ately find that 


(X , (n) = — 
i V 




n 

1=1 

Jl^^n 


'’■n ' ’‘d 


i ••= 0,1,2, •••, V - 1 


which completes the derivation. 
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Introduction 

Section VIII of the above referenced report introduced the reader to the 
constraim d moment problem wherein a solution to the classical one variable 
moment problem is sought subject to the constraint that some of the probability 
mass points are fixed a priori. While the constrained moment problem was posed 
in its general form (see Eqs. (8-2) - (8.4)), only the solutions for a few 
special cases were actually discussed. These special cases included the situa- 
tions where either one or both of the end mass points of the approximating 
probability density function (pdf) were fixed. 

Often one is interested in cases where it is desirable to fix, a priori, 
one or more of the interior mass points of the approximating pdf. (Examples 
where this situation is applicable will be discussed in the next section.) The 
solution to this more general problem was not discussed in the original report 
and is the subject of this addendum. To avoid unnecessary duplication, it will 
be assumed that the reader is familiar with the material in Section VIII and 
thus r.i erence to key equations in that section will be made wherever convenient. 
As such, the material .’iscussed here should be considered as if it was originally 
integrated into the report with the only reason for not doing so being that it 
was not available at the time the report was issued. 

The General Constrained Moment Problem 

Recall that the motivation for solving the general unconstrained moment 
problem was the evaluation of 



f(x) p(x) dx 


( 1 ) 


1 





where f (x) was "arbitrary" and p(x) was known only in terms of its first N+1 
moments 

00 

= Ejx^l = J p(x) dx; k = 0,1,2,..., N (2) 

~00 

Although never explicitly stated , f(x) was assumed to have no jump discontinuities 
since otherwise the approximate evaluation of (1), namely, 

V 

E|f(x)} = ^ 

1=1 

where the mass points x ; v and probability weights co. ; v 

are determined from the unconstrained solution of the moment problem, would not 
yield the most accurate solution. Rather, what would be desired in this situa- 
tion would be an approximating solution of the form 

V P 

6{fWl - Y, 2 

1=1 j=l 

where y^j y^? y at*e a set of fixed points corresponding to the locations 

12 p 

of the p jump discontinuities in f (x) . The solutl-"- jo this problem is clearly 
an application of the general constrained moment problem described by Eqs, (8.2) 

- (8.4) of the referenced report. 

Before proceeding to the solution of this problem, we cite a simple 
example of where an approximating evaluation such as (4) might be of use. Con- 
sider the problem of evaluating the amount of probability P in a given closed 
interval [a,b] of the pdf p(x) which is knowm to exist over the doubly infinite 


2 . 

J 


) + f(Yj ) 


(4) 


2 


interval but whose form is known only in terms of its N+1 moments as in (2) • 
Thus, we wish to evaluate 


P 


p (x) dx 


which can be written in the alternate form 


( 5 ) 



f (x) p (x) dx 


( 6 ) 


where 


f(x) 


1 ; c < X < b 


— ; x=a, x-b 


0 ; otherwise 


(7) 


Using (4) , the approximate evaluation of (6) would have the form 


A 

P 


2 “l + "2 

i=l 


( 8 ) 


where we have employed the constraints finding the solution. Note 

also that ^ ^ corresponds to the dimension of the set of unconstrained points 
which fall in the open interval (a,b) , 

With the above as motivation, we now proceed to discuss the solution to the 
general constrained moment problem. 

Let us start as before by considering the special case of p=l, where, how- 
ever, the unconstrained mass points x^, x^, are net necessarily all 


3 



required to lie above or below the constrained mass point Thus, our goal is 

to find the fewest points , x^y and probabilities 


Pr(X ^ x^) ^ J i = 1,2, V 

Pr(X = y^) = 2 ^ 


(9) 


that yields the given moments 


V 

; k = 0,1,2 N 

i=l 



( 10 ) 


Let q^, q^, and q^ be real numbers and define the polynomial 


q(x) = X + X 


(11) 


Next consider 




‘’l\+l 




Jl=l 






+ z 










k 

Z I ^ 


(=‘t) 


+ z y 

r 1 


£=1 


'K) 


We now require that q(x) of (11) satisfy the conditions 

q(yj^) = 0 


and 


q(x) > 0 for all 


( 12 ) 


(13) 


(U) 


4 


Then, using (13) and (14) in (12) gives 


V >0 


■>0^ “iVl + ’ 2 V 2 ■ 2 ’■« 


Next note that for k=0, (15) becomes 


“o'‘o ’ 2>‘2 ■ S “2"(’ 


and thus, dividing (15) by (16) produces 


Vk *'2’‘k+l 

Vo + “l"l + '>2''2 


“l-^N \ k 

V 1 a 


SL=1 V' 

I / (D.q/x \ 

\ ^ J ' j) 

\j=l 


Defining 


* A 'll’^k+l + ^2\+2 

' V0’^^1^'*'V2 


u)„ ^ — i > 0 

£ V 


j=i 


5 


we arrive at the new unconstrained moment problem given by 

V 

P./ ^ .... N-2 (20) 

£=1 

to which wo can apply our usual solution (Section III of the referenced report) 
to find x^, x^, and o)^ , , . . • , . Once having solved this uncon- 

strained moment problem, we can obtain our desired results, namely (9), from 
(16) and (19) as 



the latter result representing the normalization condition as in (8.4) of the 
ref erenced report . 

Let us now examine some special cases when X is a random variable bounded 
between a and b. 



which clearly satisfies (13) and (14). 


( 22 ) 
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Case II: = b (Constrained upper end point) 

Here we choose 


q(x) = b = X (23) 

which again clearly satisfies (13) and (14). These two cases are identical to 
the corresponding two cases given as examples in Section VIII of the referenced 
report. 

Case III: a < y, < b 

•^1 

The appropriate choice for q(x) is now 

q(x) = (x - y^j^ = - 2y^x + x^ (24) 

i,e,, a double root at x = y^. Comparing (24) with (11), we can immediately 
identify that 


■>0 ' 




(25) 


^2 = ^ 


Finally, substituting (25) into (18) and (21) gives the specific desired results 



(26) 
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(27) 


The previous results can easily be generalized to the case of two or more 
point constraints. Specifically, we are now trying to solve the most general 
problem described by Eqs. (8.2) - (8.4) of the referenced report where the p 
constrained mass points may or may not include the end points. 

To solve this most general case define the polynomials 

qj (x) ; j = l,2j . . . , p (28) 


where q.(x) is the smallest degree polynomial that satisfies 



(x) > 0 


for all X 7 ^ y. 


(29) 


Note that if y^ is an interior point then (x) will be second degree, whereas 
if y. is an end point, q.(x) will be first degree. Next, define 


P 

Q (x) = 9 j (x) 

j=l 


= Qr 


+ Q^X 


+ Q X 
m 


m 


(30) 


I 


I 


8 




where tn is the sum of the degrees of the polynomials q^Cx), q 2 (x), .... qp(x) 
Analogous to (12), consider now 


i=0 1 «.=! 






But from (29) and (30) , we have 


Q(y^) =0 ; j = 1,2, . . . , p 


Q(x^) '-0 ; il = 1,2 V 


Hence, (31) simplifies to 




V >0 

“k+i ■ E *11^ 


Evaluating (33) at k.=0, and dividing (33) by this result gives a relation 
analogous to (17), namely. 




k+i V 


Ev. "ME “3 


V i Jl 


.Q(x.) 
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Again defining the new moments 


and weights 



( 35 ) 


(36) 


we arrive at the unconstrained moment problem given by (20) where the largest 

value of k is now N-m. Once having solved this unconstrained moment problem 

for X,, x_, X and oj, , , . . . , (d , we can obtain our desired results 

12 V 12 ’ V 

from (33) with k=0 and (36) as 



( V ^ 



nQ--. 

^ J 2 J 



\j=0 / 


■"Jl Q(x^) 


(37) 


and z_ , z_, z which are the solutions to the set of linear equations 

12 p 



(38) 


Let us again examine some special cases when X is a random variable 
bounded between a and b* 
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o 
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Case IV: P = = a, « b (Constrained End Points) 

For this case, we choose 


Q(x) = (x - a) (b - x) = -ab + (a + b) x - (39) 

This example is identical to Case III in Section VIII of the referenced report. 

Case V: p = 3 , = a, a < < b == b 

Here we choose 

Q(x) - (x - a) (x - y^)^ (b - x) 

- -aby^^ + (2aby^ + (a -h b) y^^j x 

- |ab + 2(a + b) y^ + y^^) x^ -f (a + b + 2y^) x^ - x^ (40) 

Comparing (40) with (30) , we may immediately identify the coefficients Q^; 

1=0, 1,2, ..., 4 and proceed to find the desired solution from (37) and (38). The 
details surrounding other special cases of further complexity are left as an 
exercise for the reader. We do, however, point out that the recursive method for 
finding the coefficients of a polynomial generated by a product of first degree 
factors discussed in Appendix A of the referenced report is particularly helpful 
in finding the coefficients of Q(x) in (30), Note that the method used in 
Appendix A to arrive at (A-7) does not require that all the factors correspond 
to distinct roots. Thus, each second degree polynomial just be 

looked upon as a product of two identical first degree polynomials. 
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